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ABSTRACT 

o 

(•~^ ■ We simulate the evolution of an initially weak magnetic field in forced turbulence for 

(N . a range of Prandtl numbers. The field grows exponentially with the Kulsrud- Anderson 

^ . k^/'^ spectrum until the magnetic energy approaches the viscous-scale kinetic energy, 

where the magnetic forces then backreact on the velocity. Further growth proceeds more 
slowly until a saturated state is reached where the magnetic and kinetic energies are 
equal, and where the magnetic energy exists primarily at the resistive scale. We discuss 
the structure of this turbulence and the extrapolation of the results to astrophysically- 



> 

QQ ' large Prandtl numbers. 

o ■ 
o 

1. Introduction 

\ Microgauss magnetic fields are observed in spiral galaxies and between galaxies in clusters 

Ph; [Zweibel & Heiles (1997), Kronberg (1994), Vallee (1998), Beck et. al. (1996)]. In cluster plasmas 

O ' the fields have coherence lengths of up to lO/cpc (Taylor et. al., 1999). In galaxies the magnetic 

' fields are ordered over the whole galaxy and have energy densities comparable to the turbulent 

• energy density. It is not known if the fields originated before during or after galaxy formation. 

Most current research centers around dynamo theories where turbulent motions amplify the field 
^ , from an initial weak seed field to its present strength and structure [Ruzmaikin et. al. (1988), 

^ I Kulsrud (1999), Zweibel &: Heiles (1997), Beck et. al. (1996)] Indeed some kind of dynamo seems 

to be the most plausible explanation of the observations. However, despite considerable progress 
over forty years the dynamo theory is not complete and thus the history of galactic and extra 
galactic fields is uncertain. The galactic dynamo (if it exists) differs from the more familiar solar 
dynamo [Mestel, Weiss] and geodynamo [Glatzmaier &: Roberts 1995] in two key ways. First the 
disc geometry of the galaxy clearly affects the magnetic field dynamics. Secondly the ratio of 
viscosity to resistivity, the magnetic Prandtl number (Pr), is of order 10^^ in the warm partially 
ionized interstellar medium but only 1 — 10~^ in the solar convection zone. We shall also discuss 
here the possibility of dynamos in fully ionized plamas such as protogalaxies and early galaxies. In 
this paper we show that high Prandtl number dynamos are profoundly different from low Prandtl 
number dynamos. 
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The galactic spatial scales and their associated timescales suggest a possible scenario for field 
growth. Elements of this scenario are suggested by the work of Kulsrud and Anderson (1992) and, 
particularly for Stage 3, Field, Blackman and Chou (1999). We take this scenario as a framework 
for discussion and not as proven. Indeed our calculations already expose flaws in the scenario. First 
consider the important space and time scales from the largest to the smallest. The galaxy itself 
is about 10^° years old, 10 — 15kpc{= \g) across and rotates with velocity ~ 200km/ s once every 
2 X 10^ years. Supernovae produce random velocities of order 10km/ s on a scale 100pc{= A/) and 
with timescales of order lO''' years. Approximately 5 — 15% of the kinetic energy at the supernova 
scale is helical [Moffatt 1978]. Without a magnetic field the kinetic energy cascades to small viscous 
scales, where Xi, ~ 0.1 — O.Olpc. The viscous eddies turnover on a timescale ~ 10^ years. 
Finally the smallest scale is the resistive scale, which is typically lO^km - truly negligibly small. 
We will assume that the initial field is very weak. Such fields can be made in many different ways, 
see for instance Gncdin &; Ferrara, & Zweibel (2000) who find that fields of order 10~^*^ Gauss can 
be made in shocks during the re-ionization phase of structure formation. We divide the growth 
into three stages. 

Stage 1. The Kinematic Small Scale Field Dynamo. During this stage the field is 
too weak to affect the velocity at any scale. Since the eddies at the viscous scale turnover the 
fastest they amplify the field first. Therefore the growth time is approximately 10^ years in our 
galaxy. This growth was first predicted by Batchelor (1950), and later Kazantsev (1968), who 
also investigated the structure of the field, and still more recently Kulsrud and Anderson gave a 
spectral theory of the process (1992). In the 1990s the kinematic theory was considerably refined 
[Schekochihin, et. al. (2001), Chertkov ct. al. (1999)] and it is well understood. Much of what is 
understood comes from the short correlation time kinematic model in which the velocity correlation 
time is assumed infinitesimal. The important features of the magnetic field evolution in this model 
are: first, the magnetic spectrum Ei,(k, t) rises as k^^^ at all k much less than the peak, kp. Second, 
for k < kp, Ej,{k,t) grows as exp (S/Ajt) at fixed k with 7 roughly the turnover rate of the viscous 
scale eddies. The peak wavenumber, kp increases as kp oc exp (3/57^) until it reaches, and remains 
at, the resistive scale (k^ l/{\jj) ^ a/t/^ with rj the resistivity). Finally the magnetic field is in 
a folded state [Schekochihin et. al., 2001], where the variation of B along itself (< |B • VBp >) 
is much smaller than the variation of B across itself (< |B x (V x B)|^ >). All these features of 
the kinematic phase evolution are seen, at least qualitatively, in the early stages of our simulations. 
Clearly the field predicted by the small scale field dynamo is on scales much smaller than the 
observed galactic field. Therefore Stage 1 is necessarily a transitory stage. When the magnetic field 
energy becomes of order the energy in the viscous scale eddies the kinematic stage (Stage 1) ends. 
In the galaxy this corresponds to a magnetic field strength of 0.1/iG. 

Stage 2. Approach to Equipartition. While magnetic forces at the end of Stage 1 can 
change the viscous scale flows they are not strong enough to affect the more energetic larger scale 
motions. The strain of these larger scale motions continues to amplify the field. One would expect 
that they act on the field in a way similar to the viscous scale eddies in Stage 1, but on a slower 
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timescale. As the field grows more of the velocity spectrum is affected by the magnetic forces (see 
Section 5.). Eventually on a timescale of a few large eddy turnover times the magnetic energy 
becomes of order the kinetic energy of the large stirring scale motions. In the galaxy the timescale 
for this is the supernova stirring time - roughly 10^ years. On this timescale the galactic rotation 
and the helical component of the turbulence have negligible effect. 

This stage is porrly understood and it is the main subject of this paper. In 1981 Meneguzzi and 
Poquet computed the turbulent amplification of magnetic fields in a ~ 1 plasma. They showed 
that without helicity the magnetic energy grew until it saturated at about 15% of the kinetic 
energy. The magnetic and kinetic spectra resemble our higher resolution, Pj. ^ 1, computations 
(for example see Fig. 3). There is no evidence from these simulations (or our own) that there is 
energy equipartition scale by scale. In the P,. ^ 1 case, the majority of the magnetic energy is in 
the scales between the viscous and (smaller) resistive scale. Almost none of the kinetic energy is 
contained in these sub-viscous scales. 

Recently there has been considerable interest in simulating Alfven wave turbulence [Maron 
& Goldreich (2001), Cho & Vishniac (2000), Biskamp & Muller (2000)]. Identical magnetic and 
kinetic power law spectra {k~°' with a = 1.5 — 1.6) are seen. These simulations start with magnetic 
energy on the large scale that is stronger or of order the kinetic energy of the forcing scale. Dynamo 
simulations must, of course, start with small fields. It is not clear that the final state is independant 
of the initial conditions - at least for times less than the resistive time of the large scales. 

At the end of Stage 2. the magnetic energy is expected to be of order the kinetic energy of the 
forcing motions. The observed field is on scales bigger than the supernova forcing scale. There is 
no evidence from our simulations that there is any large scale field at the end of Stage 2 - indeed 
the majority of the magnetic energy is in the subviscous scales. 

Stage 3. Large Scale Field Growth. 

The final stage of the magnetic field evolution is growth of the field on the largest scale, the 
galactic scale. It is believed [Poquet et. al. (1977), Meneguzzi &: Poquet (1981), Field et. al. 
(1999), Brandenburg (2000)] that the helicity of the turbulence plays a key role in this inverse 
cascade. Indeed Field, Blackman and Chou (1999) have argued that this stage is very similar to 
the kinematic {a - co) mean field dynamo [Parker (1979), Moffatt (1978), Ruzmaikin (1988)]. The 
estimated timescale for the mean field dynamo in the galaxy is 2 x 10^ years [Field et. al., (1999), 
Kulsrud (1999)]. This timescale is controlled by the slow rotation of the galaxy that gives uj and 
drives the helicity in the turbulence to give a. The fraction of the turbulent energy in the galaxy 
that is helical is about 5 — 15%?? Recent simulations by Maron and Blackman (2001) suggest that 
large scale field growth at such low helical fractions is not possible (although it is at 100% helicity 
[Brandenburg (2000)]. They also show that the fractional helicity needed to see large scale field 
growth increases with magnetic Prandtl number. There is an ongoing debate about the effect of 
small scale fields on the mean field dynamo - several authors have claimed a quenching (suppression) 
of the alpha effect [Hughes & Cattaneo (1996), Gruzinov et. al. (1996), Bhattachargee (1996)]. 
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Blackman and Field (2000) have pointed out the important role of boundary conditions that allow 
the magnetic helicity of one sign to be discarded. Without such boundary conditions they claim 
that mean field growth is on a ncgligably slow resistive time [Gilbert]. The field we obtain at 
the end of stage 2 is highly striped (i.e. it reverses in sign on a small resistive scale). Since the 
mean field dynamo motions will amplify both the positive and negative stripes one would expect 
amplification of the small scale field to dominate. However we will not consider this stage in any 
detail in this paper. It is clear that this stage is dependant on the geometry and physical conditions 
in the galaxy and therefore simulations in simplified homogeneous models may be missleading. 

There are four key results from our simulations. First, the energy of the magnetic field grows 
to a saturated level equal to the kinetic energy of the stirred fiow. Second, the magnetic field energy 
is contained in the small resistive scales. Third, the field is in the form of long thin folds (stripes) 
where the variation of B across itself is much faster than the variation of B along itself. Fourth, the 
predominant straining and folding of the field in saturation comes from the stirring scale motions. 



We investigate MHD turbulence in situations of interest to the origin of the galactic magnetic 
field. Various regimes exist depending on the state of the dynamical medium, specifically on whether 
it is charged or neutral. In an ionized plasma, if the magnetic field is sufficiently strong, charged 
particle transport is confined to magnetic fieldlines and the kinematic diffusivity has the anisotropic 
Braginskii form (section 2.1). If neutrals are present, viscosity is isotropic. If additionally the 
magnetic energy density is sufficiently high, ambipolar diffusion can arise from the differential ion- 
neutral motion. Magnetic diffusivity in all cases is due to resistivity. We also consider the artificial 
but illustrative problem involving ordinary viscosity and resistivity (section 4). We consider these 
regimes in the context of the galaxy and protogalaxy, where kinematic diffusivities greatly exceed 
magnetic diffusivities. 

The equations of magnetohydrodynamics (MHD) are 



2. Equations 



p {dtv + V • Vv) = -V P+— -V-n + — B • VB + pvV\ 




(1) 



5tB = V X (v X B) + T^V^B 
dtp + V ■ (pv) = 0, 
dtc + V ■ (cv) = fcV^c 



(2) 



(3) 
(4) 
(5) 



V-B = 



Jl = -2,uu{BB - ^I){BB - h) 



: Vv 



(6) 
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V fluid velocity 

p fluid density 

u momentum diffusivity 

c passive scalar density 

n Braginskii pressure 



B magnetic field 

P fluid pressure 

r] magnetic diffusivity 

Vc passive scalar diffusivity 

B Magnetic unit vector 



Table 2: Dynamical variables 



n is the Braginskii tensor, which applies when charged particle transport is confined to fieldlines. 
Ambipolar diffusion can optionally be added if the fluid is only partially ionized, and if magnetic 
forces are strong enough to generate significant differential ion-neutral velocities. 

We simplify equations 1-4 for applications in this paper. The magnetic field is measured in 
velocity units with b = B/\/47r. Incompressibility is assumed throughout, so we set p = 1. We 
these modifications, equations 1-4 transform to 

dt-v = -V • Vv - V(P + y)-V-n + b-Vb + (7) 

5(b = V X (v X b) + ryV^b, (8) 
V • V = 0, V • b = 0. (9) 

dtc + yr-Vc = i^cV^c. (10) 

Incompressibility is implemented by first calculating dfV without the fluid and magnetic pres- 
sure terms, and then by projecting the resulting Fourier components transverse to k. This procedure 
accounts for the divergence introduced by the Braginskii term. 



2.1. Tensor viscosity 

Particle transport in a fully ionized medium is confined to fieldlines if the collision length 
exceeds the cyclotron radius: 

1/2 3 

B ^ ^Z!?|_^T-=^/2 ^ 2.7-io-%r3/2 

This condition is easily satisfied except possibly in the early protogalaxy where the field may have 
been very weak. For ionized plasmas, we use the Braginskii pressure (Eq. 6) in place of the 
Laplacian viscosity. 
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3. Simulation 
3.1. Scales 

Wc state our assumptions for the parameters of galactic turbulence in table 4, followed by a 
discussion of our inability to directly simulate them. 



Table 4-' Galactic turbulence parameters 



Quantity 


Symbol 


Galactic 


Protogalactic 


Magnetic coherence scale 


A. 


10^2 cm 




Forcing scale 


A/ 


10^1 cm 


10^^ cm 


Forcing velocity 


^/ 


10^ cm/s 


10^ cm/s 


Viscous scale 




10^'^ cm 


10^2 cm 


Resistive scale 


\ 


10^2 cm 


10^2 cm 


Viscosity 




1020 cmVs 


1029 cm2/s 


Magnetic diffusivity 


V 


lO^'i cmVs 


10^° cm2/s 


Prandtl number 


P 


10^ 


1019 


Neutral free path 




10^^ cm 




Ion free path 


\i 


10^3 cm 


1022 cm 


Temperature 


T 


10^ K 


10^ K 


Proton density 


n 


lcm~^ 


10-2cm-3 



The range of scales involved in the galactic dynamo poses the chief obstacle to simulation. A 
spectral grid of 256^ elements delivers a scale range of ~ 50, which is enough to capture self-similar 
dynamics or to study transitions between regimes occuring at different scales. The phases of the 
dynamo can be simulated separately and then spliced together to construct a complete picture. 
Limitations exist of course, most notably when we consider the growth of a large scale field in 
the presence of a tangled small scale field. In the galactic setting, these scales differ by orders of 
magnitude. 

In the high-Prandtl linear regime, the dominant magnetic structure exists between the viscous 
and the resistive scales. We find that at least an order of magnitude of separation between them is 
necessary to capture the k^^"^ Kulsrud- Anderson spectrum. In our simulations, a Prandtl number 
of at least 2500 is required. 

Our principal concern regarding the nonlinear stage is the final state of the magnetic field in 
the limit of large inertial range and Prandtl number, and how long it takes to get there. Specifically, 
is this field dominated by large or small scale structure, or equivalently how does the spectral index 
compare to —1? Finally, does the result for larger Prandtl number differ from that for P = 1? 

In the nonlinear stage, the magnetic field is strong enough to backreact on the turbulence. 
This first occurs at the viscous scale, then at successively larger scales in the inertial range as the 
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magnetic field grows. The objective is to study the backreaction on the inertial-range structure when 

the magnetic energy is at the resistive scale. It is important that the forcing occur at a substantially 
larger scale than the viscosity so that it doesn't interfere with the magnetic backreaction. Thus we 
require « X^, « Xf. In a 256'^ spectral simulation, the effective scale range is ~ 50, and so 
we can have a factor of ~ 7 separating each scale. To study the role of viscosity and the inertial 
range, we ran 5 simulations with X^, engineered to have a sequence of values from to A/. 

3.2. Computational scales 

Table 5 defines the scales that occur in the simulation of the turbulent dynamo. When the 
Prandtl number is greater than or equal to unity, the scales have the ordering Aj > Aj^ > A,,. The 
dynamics occur between Aj and Xa, and so we set Xf ^ 1 and A,^ A^, where the simulation volume 
is a cube of side length 1 and A^ is the smallest resolved scale. We ran a sequence of simulations 
with Xi, varying between A/ and Xa- A forcing power of unity at the outer scale maintains an RMS 
outer scale velocity of ~ 1. 



Table 5: Computational scales 







Velocity at scale A 


bx 




Magnetic field at scale A 






Kinematic diffusivity 


V 




Resistive diffusivity 


A. 




Viscous scale 


Xrj 




Resistive scale 


A/ 




Forcing scale (=1) 


V/ 




Outer scale RMS velocity (=1) 


A. 




Shear scale 


Vs 




Shear velocity 


Xa 


= 3/N 


Aliasing (resolution) scale 


iV3 




Grid size 




~ A/Zv/ 


Forcing timescale 


ts ~ Xs/Vs 




Shear timescale 


x± 




Transverse folding scale 


A|| 




Longitudinal folding scale 


f-- 


= A||/A^ 


Fieldline folding factor 


s = k/{2Tr) = 


l/A 


Wavenumber 




During 


the weak magnetic field regime. 


the following 


scalinj 


gs enable us to set A,^ and A^ by 



varying ly and rj. Below the forcing scale and above the viscous scale, the velocity has the form of 
a Kolmogorov cascade with vx Vf (A/A/)^/'^ The eddy time is t ^ X/vx- The viscous scale arises 
from equating the eddy time with the viscous time. 

A. ~ {i^x'/'/vff' (11) 

Below the viscous scale, the velocity has the form of a uniform shear with the same timescale 
as the eddy time at the viscous scale: t^, ~ X^^i'^^'^vJ^^^ . Equating this with the resistive time, 

trj ~ A^/ry, we obtain the resistive scale A^ ~ i^^'^y^j'^v^l'^v^^^'^ . We engineer 77 so that A^ is equal 
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to the grid scale Xa ~ ^f/^j where N is the number of grid elements on each cube edge. This 
establishes a condition on as a function of u: 

x3/2 3/2 
A/ vJ 

Equating the viscous and resistive times, the Prandtl number is related to the ratio of scales as 

In designing a simulation with a specified Prandtl number greater than or equal to unity, one uses 
the above scalings to select u and rj. For studies of the kinematic backreaction, it is desirable to 
have Xi, small enough so that a true kinematic inertial range exists, say Aj, ~ A//8. 

Once the magnetic field is strong enough to backreact on the velocity, the shear scale is set 
by the eddies which have the same energy density as the magnetic field. In the saturated state, 
the magnetic energy grows to the forcing energy, and the shear time increases to the forcing time. 
Simultaneously, the resistive scale increases by a factor of (Xf/Xi,)^^^ from its value in the linear 
regime. 



3.3. Code 

We withhold for now all details about the code not necessary for understanding the results, 
which we put off until section B. The equations of MHD are solved spectrally and incompress- 
ibly, and with periodic boundaries. A standard 256^ grid with spatial dimension 1^ has Fourier 

wavenumbers s = k/(2iT) extending from -85 to +85, where wc have employed the 2/3 aliasing 
truncation. Wavenumbers and physical scales are related by Xk = X2its = 2tx. Viscosity and resis- 
tivity are of the fc^ type [yV'^v and ifs/'^b) unless the Braginskii viscosity is invoked. Finally, define 
the one-dimensional kinetic and magnetic energy spectra as 

E^ = j E^{s)ds Et, = J Eh{s)ds. (14) 

The parameters for each simulation are given in table 7. 



3.4. Uniform field and helicity 

The simulations in this work have zero mean magnetic field and are forced with zero mean 
kinetic helicity. Fractional fluctuations in the kinetic helicity exist at the level of 10 percent. The 
magnetic helicity is initially zero and subsequently fluctuates about zero at an amplitude of 5 
percent of the maximum potential magnetic helicity. 
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Our simulations have no mean magnetic field or helicity, to establish what large-scale can be 
created in their absence. This is the nonlinear kinematic dynamo problem. If indeed the final 
state is dominated by small-scale magnetic energy, then the necessity of other mechanisms such as 
helicity or Keplerian shear may be invoked for the next stage of the dynamo. 

3.5. Folded fieldline structure 

The magnetic field, in the course of amplification and entanglement, develops structure not 
identifiable in the power spectrum. Define 

k1 =< (b • Vb)2 > / < 6^ > (15) 
kl=<{h-V xhf > / <b'^ > (16) 

Define a fieldline folding factor / = fcx/fcy, which parameterizes the magnetic tension per 
energy and the unwinding time of fieldlines. Also define a measure of the magnetic scale based on 
the power spectrum: 

kl = J k'^Eb{k)dk/ J Eb{k)dk (17) 

We normalize with < 6^ > instead of < >^ because it has the same kurtosis statistics as 
the squared magnetic force terms. In the simulations, k±/kb ^ 0.58 with fluctuations of 3 percent, 
whereas the kurtosis varies by 200 percent (Schekochihin, et. al., 2001). Therefore, either k± or kj, 
may be used to define the magnetic scale. 



fe^ =< (b X V X b)2 > / < 6^ > 



A;^=<(V^) >/<b-> 



3.6. Timescales 

We identify timescales for resistivity (t^), vorticity {tyj), and shear (tg) by defining: 

tr, = -i- = 27r < (V X vf >~^/^ ts =< > (18) 

tr) reflects a balance between shear and resistivity, tw corresponds to the eddy shear time and 
the magnetic growth rate during the linear regime, tg is indirectly linked to the shear time through 
tg b"^ v1 ^ Xs/vg, where b is the saturation magnetic field, and As and Vg correspond to the 
scale and velocity of the dominant shear. 
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3.7. Viscous scale 

We define the viscous scale \,y to correspond to the peak of the function k^E^{k) during the 
linear regime, which reflects the dominant scale of the shear. 

3.8. Spandex waves 

MHD interactions exist which are non-local in k space. For instance, large-scale shear can 
transfer energy directly to small-scale magnetic fields. Alfven waves are an oscillatory phenomenon 
where small-scale perturbations interact with a uniform magnetic field. A similar oscillatory phe- 
nomenon exists in a tangled magnetic field. Here, a large-scale velocity perturbation generates 
a backreaction in a small-scale tangled field. Locally, magnetic forces act in all directions, but 
spatially averaged, they tend to act collectively to oppose the original large-scale perturbation. 

Define a Lagrangian displacement field C, = ae*('''^~'^*\ with k • a = and v = dtQ. Let 
b = bo + a(bo • \i.)ie^^^'^~^*\ where bo is the static non-oscillatory part with an assumed scale of 
much less than k~^. Any other time evolution in v and b is neglected, as well as viscosity and 
resistivity, b satisfies the induction equation b = / dthd = b • VC- Returing to the Navier Stokes 
equation, < dtv >= -c^V.^^'^''-'^*) =< b • Vb >=< bo • Vbo > - < (bg • k)^ > e^(k-^-^*). An 
average is taken over scale A;, which implies < bg • Vbo >~ and < (bo • k)^ >~'< (bo)^ > 
C, has an oscillatory eigenmode with phase speed < (bo)^/3 >^/^. This is equivalent to the Alfven 
speed if we only consider the averaged magnetic field component along k. 

4. Exponential growth of a weak magnetic field 

The magnetic field is initially weak during the linear stage, and the kinetic spectrum has the 
Kolmogorov form. Magnetic fields grow exponentially at the rate of the turbulent shear. Since the 
viscous-scale eddies shear the fastest, magnetic growth proceeds at the viscous timescale ti,. 

The resistive scale of the magnetic field is determined by a balance between shear growth and 
resistive decay: tg ~ We identify the resistive scale with the transverse scale Aj^ (equation 

15). Prom the beginning to the end of the nonlinear stage, A_l increases by a factor of (A//Ai,)^/^, 
as is observed in table 6. 

4.1. Kulsrud- Anderson theory for the linear regime. 

Kulsrud and Anderson (1992) (also Kazantsev 1968, Gruzinov, Cowley, Sudan 1996, Schekochi- 
hin, Boldyrev, Kulsrud 2001) found that a dynamically weak magnetic field grows exponentially 
with a fc^/^ spectrum terminating at the resistive scale. The magnetic spectrum evolves according 
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to 



dtEbik) = I {k''dlEf,{k) - 2dkEb{k) + 6Ef,ik)) - 2^Eb{k) (19) 



which has the solution 

£;,(fc,t)~fe3/2^o(^fe)e(3/4)7* (20) 
where 7 ~ / kEy(k)dk. Ey{k)k^^'^ is the shearing rate at scale k. 

We ran simulations starting from a weak magnetic field for a sequence of 5 viscosities, all with 
Prandtl number > 1. Their ID numbers are Al through A5 (table 7), and their energy and spectral 
evolutions are shown in figures 1 and 2. Magnetic energies grow exponentially in all cases until 
the nonlinear stage is reached. The growth rates tg are given in table 6, where they are expressed 
in the form Eh(t) ~ e*/*^. The A;^/^ exponent is seen only when Pr > 2500, as with simulation Al, 
where sufficient range exists between the viscous and resistive scales (figure 2). The k^^"^ spectrum 
is especially robust in the simulation with zero resistivity (simulation AO, figure 2). This simulation 
is not without concern for its physical validity, a matter which we address in sections 5.5 , 5.8 , and 
5.7. 



Table 6: Energies and timescales 



Sim 








•n 


Ey 


Ey 














s± 


s± 












linear 


sat 


sat 


linear 


linear 


sat 


linear 


sat 


linear 


sat 


Al 


5 


10-2 


2 




0.2 


.16 


.17 


1.41 


1.5 


1.55 


4.38 


8.8 


17.0 


12.0 


A2 


5 


10-3 


1 


10-4 


0.75 


.38 


.32 


.59 


.45 


.62 


1.50 


3.5 


13.0 


8.5 


A3 


3 


10-3 


1 


10-4 


0.8 


.40 


.32 


.46 


.35 


.47 


1.26 


3.5 


14.2 


8.5 


A4 


1 


10-3 


1 


10-4 


0.9 


.60 


.22 


.36 


.25 


.30 


.78 


3.1 


18.0 


9.0 


A5 


4 


10-^ 


4 


10-4 


1.5 


.60 


.16 


.65 


.12 


.23 


.30 


1.8 


14.6 


6.0 


Bl 


5 


10-2 


1 


10-5 




.16 


.17 






2.0 




12.9 




14.0 


B2 


5 


10-3 


4 


10-5 




.3 


.4 






.69 




4.4 




12.0 


B3 


3 


10-3 


4 


10-5 




.3 


.3 






.57 




4.4 




12.0 


B4 


1 


10-3 


4 


10-5 




.5 


.3 






.31 




3.5 




13.4 


B5 


4 


10-^ 


1 


10-4 




.55 


.3 






.20 




3.4 




8.6 


B6 


1 


10-^ 


1 


10-4 




.55 


.3 






.11 




3.1 




9.0 



"Linear" and "sat" denote the linear and satuated states, averaged over a suitable length of time. Ey 
and Efy are the kinetic and magnetic energies, tg is the exponential magnetic growth time during the 
linear phase, tyj is the vorticity time, trj is the resistive time, and s± is the magnetic wavenumber. 
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The 256^ simulations were run for half of a crossing time, long enough to establish the saturated 
spectra and sy, but not long enough to accurately determine E^, E^,, and tyj. 

The magnetic kurtosis, < 6^ > / < 6^ is observed to rise from ~ 3 to ~ 15 during the 
hnear regime, and then to return to ~ 3 in the saturated state. These results are discussed in 
Schekochihin, Cowley, Maron, k, Malyshkin (2001). 

For each value of the viscosity, the resistivity is assigned the smallest value such that magnetic 
energy is destroyed by resistivity rather than dealiasing. 




T)=4 lO-* 

7)=10-< 



10- 
10- 

10- 



E^(t), v=A'W 

E^(l), 1^=1-10 
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Figure 1. Left: The exponential magnetic energy growth for a sequence of viscosities. Parameters 
for the simulation (Al through A 5) are given in table 7. Right: The evolution of the magnetic 
spectrum for simulation A3. Numbers indicate times. The spectrum at the latest time is in the 
saturated state. 
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T] = 2 ■ 10~^. Right: The same quantities for simulation AO, u = 5 ■ 10"^, rj = 0. The Prandtl 



-13- 



number is undefined, however for this resolution and viscosity, any value above 10 is functionally 
equivalent. 

5. Magnetic saturation 

Saturation is the long-term asymptotic state of MHD tmbulence in the absence of a mean 
magnetic field and with zero mean kinetic helicity. We postulate that the dynamics are character- 
ized by four principles: 1) The kinetic and magnetic energies are approximately equal. 2) Magnetic 
energy forward cascades according to the forcing timescale. 3) Magnetic fieldline folding and tan- 
gling prevents the field from unwinding. 4) The magnetic field preferentially aligns with the neutral 
shear frame. These principles lead to a magnetic spectrTim of the form that is independent of 
viscosity. The spectrum is terminated at a resistive scale that is set by a balance between forcing 
and resistivity: ~ (i/ry)^/^. 

The fastest kinematic shear is due to the viscous-scale eddies if the kinetic spectrum is more 
shallow than fc^'^, otherwise it is due to the forcing-scale eddies. In the linear regime, where the 
kinetic spectrum has the form fc®/^, the viscous-scale eddies shear fastest. In this regime, shear 
functions to transfer energy from the velocity to the magnetic field, and so the magnetic field grows 
at the rate of the viscous-scale eddies. 

In the nonlinear regime, the magnetic backreaction opposes shear that is less energetic than 
the magnetic field. This effect can be nonlocal in Fourier space, where small-scale magnetic fields 
can oppose large scale shear. In the saturated state, where the kinetic and magnetic energies are 
equal, only the forcing-scale eddies are energetic enough to shear the magnetic field. Shear exists 
on smaller time and energy scales, yet we argue that some fraction of it results from a transfer of 
magnetic energy to the velocity, rather than the other way around. Or, energy can be cyclically 
exchanged between v and b, such as for spandex waves (section 3.8). 

The magnetic field has a reduced tension per energy compared to what would be expected 
from random phase structure. The reduction is typically by a factor of 10 for a 128^ simulation, as 
evaluated by the parameterization in section 3.5. Therefore, less generation of small scale velocities 
by the unwinding of small scale magnetic fields occurs. The reduction is great enough so that the 
magnetic energy cascades from the forcing scale to the resistive scale without significant loss to 
unwinding. This phenomenon is exhibited in the tension release simulations of section 5.4, and is 
reflected by the independence of the magnetic spectrum on viscosity. 

In the simulations with an initially weak magnetic field, the magnetic spectrum is dominated 
by small scale energy at the beginning of the nonlinear stage, and in subsequent evolution it stays 
this way. One may ask if the result would be different if we started instead with a large scale, 
large energy field. Any small scale field that is formed should because the fieldline structure is 
not yet folded and the tension is strong. We observe instead that the magnetic struture evolves to 
a small-scale dominated state identical to the result of a weak initial field simulation. The likely 
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reason for this is that fieldhnes tangle, and even if they are energetically able to unwind, topological 
constraints prevent them. In general, we have never observed magnetic hysteresis in any of the 
simulations. 

Figure 3 illustrates the nonlinearly saturated spectra for a sequence of viscosities, all at fixed 
resistivity. The magnetic spectra all fall in the same place, with the exception of the one having 
the highest viscosity. This suggests that the kinematics at any scale other than the forcing scale 
are irrelevant to the cascade. The kinetic spectra also have greatly reduced magnitudes compared 
to the magnetic spectra at high k. 
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Figure 3: The saturated spectra for a sequence of viscosities and fixed resistivity. Lines with points 
indicate the kinetic spectra, and lines without points indicate the magnetic spectra. The simulations 
in the left and right panels have 128'^ and 256^ resolution, respectively. The ID numbers for the 
128^ simulations are Al through A5, and for the 256'^ simulations they are B2 through B6. 



5.1. Magnetic growth during the nonlinear phase 

During the linear stage, the magnetic energy is less than the viscous-scale energy (5 < Vy), and 
the velocity spectrum has the Kolmogorov form. The fastest shear {tg) is due to the viscous-scale 
eddies: ~ Xi,, Vg ~ Vi,, and tg ~ ti,. The nonlinear stage begins when h ~ Vi,. Thereafter, growth 
slows as the shear velocity grows according to Vg ~ h. At and above the shear scale, the velocity 
still has the Kolmogorov form, and the shear time is tg ^ Xg/vg ^ XfVg/v^ ^ Xfb^/v^. Growth 
ends when b ^ Vf and tg ^ tf. 

The resistive scale of the magnetic field is determined by a balance between shear growth and 
resistive decay: tg X^/t]. We define = A^ (equation 15). From the beginning to the end 
of the nonlinear stage, Aj^ increases by a factor of (Xf / X^)^^^ . The data in table Gsupports this 
scaling, although comparison is inexact because the forcing velocity, and hence the shear time. 
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change slightly from the linear to the saturated state. 

The magnetic enegy grows during the nonlinear stage in our 3 dimensional simulations, whereas 
in the 2 dimensional simulations of Kinney et. al. 2000, it decays. In 2D, they found that an initially 
weak field grows exponentially at the resistive scale until nonlinearity occurs, after which the field 
decays. 

5.2. Timescales and the neutralization of shear 

According to the backreaction hypothesis, shear motions below the forcing scale have a dimished 
role in cascading the magnetic field. If true, then the vorticity time reflects some shear motions 
not associated with the cascade of the magnetic field, whereas the resistive time more accurately 
reflects the magnetic cascade time. We therefore expect a greater change in from the linear to 
the saturated state than for tyj, and we also expect to be more approximately constant as a 
function of viscosity than tw Both of these conclusions are well supported by the data in table 
6. is not expected to be exactly constant in the saturated state because the larger viscosities 
influence the forcing-scale velocities, and hence Ey and the shear rate. 

5.3. Large initial magnetic field 

Kulsrud- Anderson predict that an initially weak field grows exponentially with a k^^"^ spectrum, 
which simulation confirms. The A;^/^ spectrum is independent of the initial spectral shape. The 
magnetic energy is therefore dominated by small-scale structure at the beginning of the nonlinear 
phase, and remains so through subsequent evolution to the saturated state. One may ask if the 
result would be different if the magnetic field were initially strong and organized at low k. We will 
find that there is no hysteresis in the saturated state from the results in this section, and in section 
5.7 on resolution refinement. 

Simulation L3 has an initial magnetic energy of unity, which is confined to modes with k/2iT < 
4. The viscosity is 3 • 10"'^, and the resistivity is 10~^, both of which are identical to simulation 
A3. The viscous scale is smaller than the initial magnetic scale. Subsequent evolution restores the 
magnetic spectrum to the saturated state of simulation A3, as show in figure 4. 

A more stringent test is to initialize a strong large-scale magnetic field which has no topological 
linkage. Such a field can resist bending by the turbulence and possibly oppose the formation of 
small scale field. If the field is topologically linked, then even though tension forces exist, they may 
be unable to release the tension. The question is, do linkages form in a strong magnetic field? 

The field of 62 = 2sin(27ra;i) is a simple example having zero mean and an energy density 
of unity, and constitutes the initial conditions of simulation L4. The initial velocity is zero, and 
the forcing power is unity. The viscosity is 10~^ and the resistivity is 10~^, which are identical to 
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simulation A4. The magnetic field evolves to the saturated state A4 after 5 time units. 
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Figure 4- Left: this figure follows the evolution of an initially large magnetic field at low k (simula- 
tion L3), until it reaches a state where the magnetic energy is dominantly at high k. The dashed line 
is the time-averaged saturated state of simulation A3, which has the same viscosity and resistivity. 
Right: This figure follows simulation L4, which is similar to L3 except that the initial field is at 
even lower k. The dashed line is the time-averaged saturated state of simulation A4, which has the 
same viscosity and restivity. 



5.4. Interaction between the velocity and the magnetic field 

The reduced magnetic tension of the satTiratcd state was previously diagnosed by noting that 
^\\/^± >> 1- A more direct test can be constructed by observing how fast magnetic tension 
generates kinetic energy. This is done by artificially setting the velocity to zero and observing the 
release of tension in the subsequent evolution. Specifically, we note how much new kinetic energy 
is generated from fieldline unwinding, and where in Fourier space it appears. Folded fieldlines 
have less tension per energy than structureless fieldlines, and hence unwind more slowly, generating 
less velocity. The random-phase magnetic field defined in appendix C serves as a reference for 
structureless fieldlines. A random-phase field is generated from a structured field by randomizing 
the Fourier component orientations while preserving the power spectrum. The structured magnetic 
field generates more velocity than the random-phase field, verifying that it has reduced tension per 
energy (figure 5). 

We draw from the nonlinearly saturated state of simulation A4 at t = 12 to initialize two 
test simulations. In simulation U4, we erased the velocity and restarted with the same viscosity. 
The resistivity was reduced to 10~^ to remove the effect of resistive magnetic energy loss. In 
simulation U4r we additionally randomize the phases of the magnetic Fourier modes. A random 
phase magnetic field has a folding factor of unity, and serves as a reference for assessing structure. 
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By comparing the kinetic spectra generated in simulations U4 and U4r, we can evaluate the degree 
of folding present in U4. 

We also observe that small scale magnetic fields generates large scale velocities, while the 
random-phase field generates small scale velocities. Specifically, the original magnetic field at 
s ~ 15 generates velocities at s ~ 3, while the random phase field generates them at s ~ 15. 
We infer that for the dynamics of the saturated state, only large scale shear is responsible for the 
magnetic cascade, even though the kinetic spectrum is more shallow than k~^. Therefore, the 
cascade proceeds according to one timescale, the forcing timescale. This also supports our claim 
that magnetic forces oppose shear that is less energetic than the magnetic field. Finally, the absence 
of production of kinetic energy beyond s = 6 establishes that unwinding of the small scale field 
is unimportant. The magnetic energy cascades from the forcing scale to the viscous scale without 
loss to the velocity. 




1 10 100 1 10 100 



k/2n k/2n 

Figure 5. Left: Simulation U4- Right: Simulation U4r. Numbers indicate times. These simulations 
are discussed in section 5.4- At t = Q, we show the kinetic spectra from before the velocity was 
artificially set to zero. 



5.5. Folding scales 

We evaluated the longitudinal and transverse folding scales Ay and A_l of the saturated states 
for a range u and rj (figure 6). For a 128^ simulation, and for rj larger than 10^^, wc find that 
A_L ~ ry^/^, which would be expected if the shear timescale is constant and shear balances 

resistivity. Ay is approximately constant and equal to the outer scale of unity. For rj < 10~^, 
decreasing the resistivity does not further decrease A_l, but it does decrease in Ay and the folding 
factor. This occurs because for these low values of 77, dealiasing of the magnetic field disrupts folded 
structure (section 5.7). We interpret the lowest attainable value of A_l as the resolution limit. 
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For a 256^ simulation, the resolution limit is reached at 77 ~ 4 • 10~^. At this resolution, Aj^ is 
smaller than for a 128^ simulation, while Ay is unchanged. Further reduction of r) does not further 
reduce A_l. 

The saturated spectra for varying resistivity at fixed viscosity are generated by a process of 
resaturation, whereby we take a saturated state, change the resistivity, and continue evolution 
until the system has reached a new equilibrium state. This process is described in section 5.7. The 
initial conditions are always taken from one of the simulations in the series Al through A5. 256^ 
simulations are initialized by doubling the grid of the appropriate 128^ simulation, resetting the 
resistivity, and then evolving to the new saturated state. 
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Figure 6: The upper and lower rows of points correspond to the longitudinal and transverse folding 
scales defined in section 3.5. These numbers come from the saturated states of simulations for a 
range of values in v and rj. Small symbols represent 128^ simulations, and large symbols represent 
256^ simulations. The square points in the figure represent a 64^ simulation with a Braginskii 
viscosity of vt = -003, and with 77 = .0004. This establishes that fieldline structure is folded in this 
simulation. 



5.6. Unwinding 

For scales large enough to be uninfluenced by viscosity, magnetic fieldlines release tension 
in one Alfven time: tA ~ Below the viscous scale, the viscous and the magnetic tension 

terms dominate in the Navicr Stokes equation: vv,u/\\ ~ 6^/A||, where Vu is the velocity at which 
magnetic structures zunwind. The viscous unwinding time is then t^ ~ \\/vu ~ We 
have accounted for the possibility of folded magnetic structure by distinguishing between A_l and 
A||. 

Magnetic energy is lost Alfvenically if Ia > tu, otherwise it is lost viscously or resistively. 
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whichever is faster. The transition between Alfven and viscous unwinding occurs at a scale Xa 
such that tA ~ tu- Assume a spectrum of Ei){k) ~ k~^, or & ~ A*^. We consider two limiting cases 
for the folding structure. If Ay ^ A/, then A^ ~ A^^A^^^. If Ay ~ A_l, then A^ ~ A^^^A^ 

At small scales, where viscosity mediates unwinding, we ask which is faster, shear or unwinding. 
If ts < tu, then the magnetic energy cascades without loss to the resistive scale, which is determined 
by a balance between shear and resistivity: A^ ~ {tsr])^^'^- The magnetic energy would then grow 
to be equal to the forcing energy. If tg > tu, then the magnetic energy must adjust itself to restore 
is ~ i« ~ We define the resistive scale to be equal to the transverse magnetic scale: A^ = A_l. 

The ratio of the unwinding to the shear time is 



ts A2 62 r? ^^^^ 

Simulation suggests that the magnetic field and the shear eddies have the same energy: 6^ ~ Ug, 
and that the folding scale is such that Ay ~ A/. We assume that the Kolmogorov scaling applies 
above the shear scale: v^Xf ~ vj-Xg. 

If tu > ts, then b ^ Vf and Xs = Xj, and consistency with (21) requires that 

A.;:|>1. (22) 

If A < 1, then unwinding balances shear, and decreases to bring about i„ ~ ~ i^. The 
shear time and the magnetic scale both contract until A^ ~ (ryi/)^/^ A^/^. The magnetic energy is 
then 6^ ~ vjA^/^. 



5.7. Resolution refinement 

We found in section 5.3 that the saturated state is free from hysteresis. There, we compared 
the final state of simulation A4, which started from weak magnetic field conditions at high k, with 
that of simulation L4, which started from strong magnetic field conditions at low k. The viscosity 
and resistivity are the same in each of these simulations. In this section, we consider the effect of 
changing the resistivity. 

In a process which we call resaturation, we take the saturated state from a finished simulation, 
change the resistivity, and then continue the simulation until a new saturated state is reached. Sim- 
ulations A4s and A5s were thus obtained from simulations A4 and A5, respectively. In simulations 
A4s and A5s, the resistivity was reduced to 10~^. 

Simulations A4w and A5w have the same u and r) as A4s and A5s, but they were started 
instead from a weak initial magnetic field with an energy of 10~^. The saturated spectra for A4s 
and A4w are identical, as are the ones for A5s and A5w. 
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The new saturated states for simulations A4s and A5s were reached in a time equal to one third 

of the forcing time. Whereas the technique of resaturation has now been validated, we can thus 
obtain the saturated states for any combination of v and r/ with little computational effort. This 
procedure was exploited to obtain the data used in section 5.5, and also to obtain the saturated 
spectra of the 256^ simulations Bl through B6. 

Simulation A4s develops only slightly more small-scale magnetic field than simulation A4. 
This is because the resistivity is small enough so that the spectral aliasing procedure contributes to 
magnetic energy loss, which destroys folded fieldline structure and increases the magnetic tension 
per energy. We show this with the tension spectrum in figure 7. The extra magnetic forces 
generate velocities which dissipate viscously, removing magnetic energy. Viscosity takes over the 
role of resistivity, diminishing the effect of lowering the resistivity in simulation A4s. 

Although not obvious in the spectra, this effect can be seen by plotting the viscous and resistive 
dissipation as in figure 8. The aliasing scale is reduced by a factor of 2 in a 256^ simulation, where 
it then does not trigger viscous dissipation. 
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Figure 7: We plot the spectrum of the tension term, b • Vb, to show the effect of lowering the 
resistivity enough so that aliasing destroys folding structure and enhances tension. The simulations 
are A4 and A^w, which both have v = 10~^. A4 has a resistivity of 10~^, and A^w has a resistivity 
of 10~^. The random phase spectra are obtained by the random phase transformation defined in 
section C 



5.8. Dissipative power 

We compare the spectral distribution of energy and dissipation for the saturated state by 
plotting each quantity logarithmically distributed in k. For the kinetic and magnetic spectra, we 
plot ln{W)kE{k), and for the viscous and resistive power, 2uhi{W)k^E^ik) and 2r? In(10)fc3£;^(fc). 
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When plotted on a linear ordinate, these curves have the interpretation that the total quantity is 
equal to the area under the curve. 

In figure 8, it is clear for simulation A4 that the magnetic energy resides small scales, and 
that the high-k cutoff is due to resistivity. Lowering the resistivity by a factor of 3.3 in simulation 
A4t does not significantly change the magnetic energy spectrum because viscosity takes the place 
of resistivity at high k. The interpretation is that in simulation A4t, the resistivity is sufficiently 
low so that aliasing contributes to magnetic energy loss. This destroys magnetic fieldline folding at 
high k, increases the tension (figure 7), and leads to more viscous dissipation. A proper study of 
the effect of lowering the resistivity requires that we simultaneously increase the resolution. 
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Figure 8. Left: We plot the energy and dissipation per logarithmic k for the velocity and magnetic 
fields. Lines without dots denote simulation A4, with u = 10~^ and r) = 10~^. Lines with dots 
denote simulation A4t, which has the same viscosity and has 77 = 3 • 10~^. Right: The energy and 
dissipation per logarithmic k for the 256^ simulation B4, which has u = 10"^ audi] = A-IQ-^ . The 
magnetic spectrum of simulation B4 farther toward the right in k than for simulation A4- 



5.9. Shear and magnetic alignment 

The symmetric shear matrix S = {diVj + djVi)/2 can be diagonalized to yield the shear 
eigenvalues {A+,Aj,A_} and eigenvectors {s+,so,s_}, where we list the eigenvalues from largest 
to smallest. Incompressibility ensures that A+ + A/ + A_ = 0. The magnetic energy growth is 
G = .5 < dtlP' >=< bibjSij > . With the magnetic field projected into the local shear eigenvector 
frame, this is G =< b\X+ + 6qA/ + 6?.A_ > . 

In a nonmagnetized fluid, let a passive scalar c have a diffusivity ( which is substantially less 
than v. is known as the Schmidt number. In the range between X^, and A^, where A^ is 

the passive scalar diffusion scale, the passive scalar cascade rate is Cc ~ / {Xi, / v,^) , leading to a 
steady-state cascade index of Ec{k) ~ k~^. 
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Consider two extremes for the orientation of b in the shear eigenvector frame. First, if b 

hes principally along sq, fieldlines are neither being stretched nor compressed, only interchanged, 
and therefore they cascade without amplification. A profile of b in a plane perpendicular to sq 
has an appearance similar to a passive scalar. In this extreme alignment, the power spectrum is 



The other extreme is for the fieldlines to be aligned with s+, so that they are amplified as 
they cascade. If each reduction in scale is accompanied by a self-similar amplification of energy, 
the resulting power spectrum index is Ei,{k) ~ k^. The geometry of fieldlines in the saturated 
state must be bracketed by these extremes, therefore the subviscous magnetic spectrum will have 
an index between —1 and 0. This is observed in figure 3, where the index is closer to —1. The 
alignment of b in the eigenvector frame is found to be preferentially along sq, which also suggests 
an index closer to —1. 

We measure the saturated-state alignment of b in the local shear frame by evaluating G2 =< 
(bibjSij)'^ >, which is the variance of the magnetic alignment. In the time-averaged saturated state, 
G = 0, and so the mean is zero. We compare this to the result obtained after random-phasing b 
(but not v). For a structured field, G2 has one half the value as for a random-phased field. This 
establishes that in the local shear frame, b preferentially aligns with sq- 



Simulation K4 is forced at 3 < s < 4, unlike the other simulations which are forced at 1 < s < 2. 

Its purpose is to determine if magnetic energy can occupy modes larger than the forcing scale. The 
result from figure 9 is that they do not. The initial conditions are taken from the saturated state 
of simulation A4 at t=19. The viscosity and resistivity are the same as for simulation A4. 
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Figure 9: In simulation K4, magnetic energy does not occupy to a significant degree the wavenum- 
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bers smaller than the forcing wavenumhers at s = 3, 4. 



5.11. Turbulent difFusivity 

For Pr < 1, sub-resistive velocities hinder magnetic growth. We establish this with two sim- 
ulations having the same A^^ and different (figure 10). The first (SI) has Xi, = A^, and the 
second (S2) has Xu ~ 0.4A,j, or Pr = 0.4. The second has sub-resistive kinetic energy while the first 
does not. The resolution of S2 is sufficient to fully capture the dynamics. The magnetic energy 
in simulation SI grows more quickly than that for simulation S2. The interpretation is that the 
sub-resistive eddies act as a turbulent diffusivity on the magnetic field. We will see elsewhere that 
turbulent diffusivity does not apply for Pr >1. 

In figure 11, we observe that the magnetic field decays for Prandtl numbers 0.2 and 0.1. 
However, it is possible that for a simulation with these Prandtl numbers, and a larger grid, that 
the field might grow instead. 
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Figure 10. Left: simulation SI, u = W r/ = 10 Pr = 1. The magnetic field grows robustly. 
Right: simulation S2, i/ = 4 • 10~^, rj 
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Figure 11. Left: simulation S3, u = A ■ 10~^, 77 = 2 • 10~^, Pr = 0.2. Right: simulation S4, 
u = A - 10"^, 77 = 4 • 10~^, Pr = 0.1. In each of these simulations, the magnetic field decays. 



5.12. Pitfalls 



Subtle dangers threaten the validity of turbulent dynamo simulations. 

For saturated dynamics, magnetic unwinding causes kinetic motions and viscous dissipation 
at scales below Xi, (section Z). Therefore, the use of a hyperviscous operator in the Navier-Stokes 
equation of the form f„V^"v, where n > 1, is not valid. If used, it would anomalously destroy 
small-scale magnetic fields, giving the false impression that the saturated magnetic spectrum does 
not extend beyond the viscous wavenumber. 

In hydrodynamic turbulence, inertial range dynamics are insensitive to whether energy is 
removed at the inner scale by physical diffusivity or by dealiasing. In MHD turbulence, dealiasing 
destroys folded fieldline structure, increasing the unwinding and viscous dissipation by a factor of 
(A||/A_l)^, which can be as large as 200 in our simulations. Therefore, the resistivity must be large 
enough to prevent magnetic energy from reaching the aliasing scale. 

The standard a — Q helical dynamo theory assumes that locally, the uniform magnetic field 
component outweighs the fluctuating component. It is of interest to know if the turbulent dynamo 
yields a magnetic configuration in accord with this assumption. We therefore emphasize zero 
helicity and high Prandtl number in our simulations. The astrophysical context has an inertial 
range of ~ 4 decades and a Prandtl number of ~ 10^^. 

The magnetic field backreacts on the turbulence as it increases in energy after the linear regime, 
slowing down the shear timescale and increasing the resistive scale according to tg ^ X^/rj. Although 
the magnetic scale increases, this should not be interpreted as the formation of a large-scale field, 
as the field is still at the resistive scale. Large resolution and Prandtl number are required to yield 
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a saturated state where the magnetic energy scale is well separated from the kinetic energy scale. 

The a — Q helical dynamo theory invokes a concept of turbulent diffusivity acting on the 
magnetic field. In this picture, turbulence with RMS velocity v at scale A acts diffusively on the 

magnetic field with an equivalent viscosity of v\. This is at odds with our results, which suggest 
instead that the magnetic field is cascaded non-diffusively, and only by velocities which have more 
energy than the magnetic field. 



6. Braginskii viscosity 



Simulations with Braginskii viscosity are qualitatively similar those with Laplacian viscosity. 
In figure 12, we see that the field grows quickly in the linear regime, and then reaches a saturated 
state that is essentially the same as that for the Laplacian viscosity simulations. The difference in 
the dynamics is that the Braginskii viscosity does not oppose fieldline unwinding, and so the Alfven 
scaling would apply for the unwinding time. We speculate that topological constraints inhibit the 
unwinding of small-scale magnetic fields. This is similar to our interpretation in section 5.3 for 
why a large energy, large-scale magnetic field ultimately evolves to a small-scale field. 
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Figure 12. Left: The evolution from the linear to the saturated state for the Braginsky viscosity 
simulation T3, with u = 3 - 10~^ and 77 = 10"^. Right: The saturated magnetic spectra for two Bra- 
ginskii viscosity simulations (T2 and T3), shown together with the Laplacian viscosity simulations 
from figure 3. The resistivity is 10~^ in each. The magnetic spectra are nearly the same for each 
type of viscosity. 



7. Acknowledgements 



We wish to thank James McWilliams, Alex Schekochihin, Benjamin Chandran, Eric Blackman, 
Ellen Zweibel, Yoram Lithwick, and Axel Brandenburg for helpful discussions. We benefitted from 



-26- 



the supercomputers at the Caltech Center for Advanced Computing Resources (10000 CPU hours) 
and at the National Center for Supercomputing Apphcations at UIUC (10000 CPU hours), and 
from their very helpful staff. 



A. Index of simulations 

Table 7: Index of simulations 



ID 


Grid 




u 






'/ 




Pr 


Notes 


AO 


128^ 


5 


10 


-2 













Al 


128^ 


5 


10 


-2 


2 


10- 


-5 


2500 




A2 


128^ 


5 


10 


-3 


1 


10- 


-4 


50 




A3 


128^ 


3 


10 


-3 


1 


10- 


-4 


30 




A4 


128^ 


1 


10 


-3 


1 


10- 


-4 


10 




A5 


128^ 


4 


10 


-4 


4 


10- 


-4 


1 




Bl 


256^ 


5 


10 


-2 


1 


10- 


-5 


2500 




B2 


256^ 


5 


10 


-3 


4 


10- 


-5 


125 




B3 


256^ 


3 


10 


-3 


4 


10- 


-5 


75 




B4 


256^ 


1 


10 


-3 


4 


10- 


-5 


25 




B5 


2563 


4 


10 


-4 


1 


10 


-4 


4 




B6 


2563 


1 


10- 


-4 


1 


10- 


-4 


1 




L3 


128^ 


3 


10 


-3 


1 


10- 


-4 


1 


Init b at fc = 4 with large energy. 


L4 


1283 


1 


10 


-4 


1 


10 


-4 


1 


Init b at fc = 1 with large energy. 


T2 


643 


5 


10 


-3 


1 


10- 


-4 


50 


Braginskii viscosity 


T3 


643 


3 


10- 


-3 


1 


10- 


-4 


30 


Braginskii viscosity 


T3b 


643 


3 


10 


-3 


4 


10- 


-4 


7.5 


Braginskii viscosity 


U4 


1283 


1 • 


10- 


3* 


1 


10- 


-5 


100 


Erase v from A4 and continue. 


U4r 


1283 


1 ■ 


10" 


3* 


1 


10- 


-5 


100 


Like U4, & with random-phased b 


K4 


1283 


1 


10 


-:S 


1 


KJ- 


- 1 


100 


F()r(:(^ at, s=3 & 4. 


A4w 


1283 


1 


10 


-,'S 


1 


10- 


-6 


1000 


Start from weak field 


A5w 


1283 


4 


10 


-4 


1 


10- 


-6 


1000 


Start from weak field 


A4s 


1283 


1 


10 


-3 


1 


10- 


-6 


1000 


Start from saturated state of A4 


A5s 


1283 


4 


10 


-4 


1 


10- 


-6 


1000 


Start from saturated state of A5 


SI 


1283 


1 


10 


-3 


1 


10- 


-3 


1 




S2 


1283 


4 


10 


-4 


1 


10- 


-3 


0.4 




S3 


1283 


4 


10 


-4 


2 


10- 


-3 


0.2 




S4 


1283 


4 


10 


-4 


4 


10- 


-3 


0.1 





The following properties are common to all simulations: The box size is (1,1,1), and the kinetic 
energy is forced with a power of unity. Forcing occurs within a sphere of radius 2 in Fourier lattice 
space, except for simulation K4, which is forced at a radius of 3 and 4. 
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B. Code 

We evolve the incompressible MHD equations in Fourier space where they take the form 
(Lesieur 1990) 

dtVa = -ikj (^6ai3 - ^) (yv^ - bpb^ - vk'^'^Va , (Bl) 
dtba = -ikpivpb^ - bpvj - rjk'^''ba, (B2) 

kaVa = 0, kj}a = 0, (B3) 

dtc = —ikjjvpc — Uck'^'^c. (B4) 

Equations Bl, B2, and B4 constitute a system of ordinary differential equations with time 
as the dependent variable and the Fourier coefficients h^, c} as the independent variables. 
We employ a modified version of the second order Runge-Kutta algorithm (RK2) to advance the 
variables in time. 

We make one departure from RK2 and treat diffusive terms with an integrating factor. Con- 
sider an equation of the form 

dtq{k)=A-^nk^''m. (B5) 

where A comprises the non-diffusive terms. Its solution, with A constant throughout the interval 
Ai, is 



We use this expression in place of El in each stage of RK2. To lowest order in Unk'^^At, cqiiation 
B6 reduces to El. However, it has the advantage that it yields stable solutions to equation B5 with 
constant A for arbitrary values of v^k'^^/S.t whereas El yields unstable solutions for i^nk'^'^At > 2 

We employ standard dealiasing according to the 2/3 rule (Canuto 1988) except with the Bra- 
ginskii and ambipolar terms, which require special treatment (section B.l). 

The code, written by Maron, is exhaustively discussed in Maron & Goldreich 2000. 

Define the one-dimensional kinetic and magnetic energy spectra as 

Ey = j E^{k)dk Eb = J Eb{k)dt (B7) 



B.l. The Braginskii term 

Bilinear terms in the MHD equations are calculated by transforming the individual fields to 
real space, carrying out the appropriate multiplications there, and then transforming the products 
back to Fourier space. This requires N1N2N3 operations using the Fast Fourier Transform (FFT) 
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algorithm; {NiN2Ns)^ operations would be needed to carry out the equivalent convolution in Fourier 
space. 

This economy comes at the price of either a 1/3 reduction in resolution or an aliasing error. 
To appreciate this, consider the ID product 



where m is any integer. The m = terms comprise the convolution, and the remainder the aliasing 
error. To avoid the aliasing error, we set all Fourier components with |s| > to zero both before 
we compute the real space fields and again after we return the bilinear terms to Fourier space. 
Truncation ensures that Fourier components of bilinear terms with m vanish. Its cost is the 
reduction of the effective spatial resolution from N to 2N/3. 

The Braginskii term contains a product of 5 fields inside the gradient, plus two divisions by 
b^. When M terms are multiplied, dealiasing should be applied for |s| < N/{M + 1), however the 
division will introduce aliased terms regardless of the location of the truncation. The appropriate 

procedure then becomes a matter of engineering. The rule we apply to decide if the computation 
is sufficiently accurate is that it should agree with the equivalent computation performed on an 
infinitely large grid. We find that a truncation of |,sl < is inadequate for the Braginskii term, 
and that a truncation of |,s| < is more than adequate (figure 13). 

To test the properties of the dealiasing truncation, we evaluate the Braginskii term for two grid 
sizes, with the larger grid regarded as the "correct" result and the smaller grid the "test" result. 
We apply the same dealiasing truncation to both grids before and after the computation before 
comparing. We then varied the location of the dealiasing truncation for different pairs of grids. In 
table 8, pairs A-B, C-D, and E-F all share the same truncation. The input fields are taken from 
simulation Z65 at t = 4. The magnetic field is dominated by small scale structure and presents a 
worst case scenario for alias error. In figure 13, we selected the Braginskii term's x component and 
plotted it as a function of x at fixed y and z for each grid pair. Test grid F, with a truncation at 
A^^/3, shows poor agreement with the larger grid E. 

We conclude that the MHD terms should be dealiased at A'^/3 and the Braginskii term at 
We implement this by computing the MHD terms on a size and the Braginskii term on a (2A^)^ 
grid. The Braginskii result is then returned to a grid and added to the MHD result. Aside 
from economizing resolution, this configuration also economizes timcstep. The timcstcp will be set 
by the grid, which is twice as large as the timestep that would be necessary for the {2N)^ grid. 



pq{s) 



_ j_ 

~ N 



= W Es' Es" p(6')5l5")e'"'^^'+^"^'/^5.,.'+."+miV, 



(B8) 



Table 8: Braginskii viscosity test simulations. 
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ID 


Grid 


Truncation 


ID 


Grid 


Truncation 


A 


64^ 


|s| < 32/6 


B 


32^ 


\s\ < 32/6 


C 


643 


|s| < 32/4 


D 


32^ 


\s\ < 32/4 


E 


64^ 


|s| < 32/3 


F 


32^ 


|s| < 32/3 



The ambipolar term is computed by first obtaining the magnetic force and then putting it in 
the place of v in the induction equation. This requres that we deahas the magnetic force before 
computing the rest of the induction term. 




-0.3 B: 32^ 



I , , I , I I I , , I 

0-5 1 

X 

Fig 13: Alias error in the Braginskii term,. Curve B is the Braginskii term computed with a 
dealiasing truncation at N/6, which agrees will with the equivalent simulation on a larger grid. The 
dealiasing truncation in curves D and F is at N/4 and N/3 respectively. The N/3 truncation leads 
to poor agreement. 

C. Random phase transformation 

Structure in the dynamical fields can be studied through comparison with their random-phase 
realizations. The random phase transformation is designed to remove the coherent structure while 
preserving the power spectrum. The transformation is accomplished by giving each Fourier mode 
vector v(A;) and b(A;) a new random orientation while preserving the magnitude, subject to the 
constraint of divergencelessness. 
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